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Emergent gravitational dynamics from multi-BEC hydrodynamics? 
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In this paper, we examine the possibility to implement some form of emergent Newtonian gravity 
in a generic multi-component Bose-Einstein condensate. Parallely to what happens for the emer- 
gence of low energy Lorentz invariance, strong requirements have to be imposed on the underlying 
condensed matter model. We will show, within a simplified model, that the presence of a global 
symmetry alleviates the problems associated to Lorentz violation, allows the presence of a long range 
potential, to which the analogue matter fields (the quasi-particles) are coupled following a weaker 
' form of equivalence principle. 
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I. INTRODUCTION 



The search for quantum gravity has lead to a variety of different approaches exploring different directions [lj , with 
different choices of foundational axioms and with different levels of success. Despite some promising candidates (like 
superstring/M-theory Q and spinfoam/LQG Q) we do not have any completely satisfactory quantum theory for the 
gravitational field. 

With respect to this quest, the celebrated results coming from the physics of black holes in general relativity [J], and, 
_j ' more recently, derivations of the Einstein equations as equation of state [5| within a fully thermodynamical treatment 
(see also [6H8[) suggest a particular perspective. In view of these results, the metric would be just a collective degree 
of freedom, pretty much like a thermodynamic variable: from this point of view, then, the direct quantization is not 
at all the right thing to do. While this idea has a certain appeal on its own, it is important to understand what are its 
limits, and ultimately, what can be said about the properties of the underlying degrees of freedom by the inspection 
of the properties of the emergent collective theory. 

Analogue models for gravity @ have pro vided some useful insights, giving us the opportunity to test with numerical 
simulations [l(| and even experiments [ill . [l2j some key features of physics in strong gravitational fields that would 
be not so easily accessible otherwise. Additionally, and perhaps most importantly, it has been possible to push a bit 
further the analysis to include the effects of some particular trans-Planckian completions on certain predictions for 
low energy physics [IH [3] and even more exotic phenomena [HI that have no direct analogue in general relativity. 

The power of the analogue models, however, cannot be of great use when we try to address fully dynamical 
phenomena: after all, fluid equations of the various sorts cannot be put into one to one correspondence with Einstein 
equations. Nonetheless, while we cannot hope to reproduce dynamical gravitational phenomena with analogue models, 
we can still hope to understand some of the symmetry principles that are at the foundations of general relativity. 

The purpose of this paper is to dig in this direction, by considering a particular condensed matter system, namely 
a Bose-Einstein Condensate (BEC). We will consider the case of a system with several species of bosons, and we will 
try to address the issue of the emergence of a sensible theory of phonons (which in the analogy will play the role of 
the matter fields), moving over the condensates (providing the notion of acoustic geometry). 

Obviously, the most crucial obstacle to the realization of a reasonable theory (within an analogue model) is the 
absence of symmetries like diffeomorphism invariance and local Lorentz invariance from the very onset. As it has 
been observed 16], the most general low energy geometrical structure describing the propagation of perturbations 
can be cast in the form of a Finsler structure, and not of a Lorentzian one. Of course, this has deep implications 
on the possibility of simulating accurately all the phenomena associated to gravity, given the role of local Lorentz 
invariance in general relativity. Nonetheless, one can hope to get some insight on some structural issues which might 
be independent from the explicit form (and content) of the equations of motion. 

In a previous work [ItJ it has already been discussed the pseudo-gravitational theory that emerges in the case of 
a single component BEC. Despite being very primitive, it is an example that explicitly shows how, from a micro- 
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scopic theory (a many-body quantum system), one can obtain some form of (acoustic) metric and a form of pseudo- 
gravitational theory in the hydrodynamic limit. Clearly, despite the obvious limitations of the system, that model 
is on the right track to implement quantitatively the otherwise vague idea that general relativity is a hydrodynamic 
theory. 

Despite not being immediately relevant for the discussion of a nonrelativistic theory, to have a form of low energy 
Lorentz invariance in this kind of models will help in bringing it closer to physically interesting situations. Therefore 
we will reconsider the emergence of Lorentz invariance at low energy in BEC, trying to extend (at least partially) the 
results of jl8|, showing that the imposition of an internal symmetry is not enough to achieve the desired result of a 
low energy Lorentz invariant effective theory, but that a tuning is still needed. However, as we will show explicitly, 
the presence of the underlying symmetry makes Lorentz violation less severe. 

For what dynamics is concerned, we will extend the analysis done in (l7| . considering some backreaction effects of 
the phonons onto the condensate. In the case of the single component BEC, it has been shown that it was possible 
to mimick some form of nonrelativistic Newtonian dynamics (with important differences that made the model totally 
useless for the purpose of realistic dynamical analogues) . Here we will reproduce the analysis in the more general case, 
with the objective of finding long range potentials, and, most importantly, a form of equivalence principle 1 governing 
the coupling of matter fields to the candidate for the "gravitational field" that we will identify. 

The plan of the paper is as follows. In section |TT] we present a Z2-symmetric model of a two-components BEC, 
where we show the emergence of a long range potential. We start the analysis of the general multicomponent case 
in section llrfl where we discuss the quasi-particle spectrum in homogeneous condensates. Section llVI deals with the 
analysis of the dynamics of small inhomogeneities in the condensates as it is induced from the underlying many body 
Hamiltonian, beyond the Gross-Pitaevski approximation. In the same section, the candidate for the analogue of the 
Newtonian potential is identified an its relationships with the inhomogeneities in the condensate elucidated. Section 
fVl is devoted to the discussion of the structure of the source terms for the potentials, while in section l*VTl we discuss 
the difficulties in obtaining a meaningful local form for the analogue of the Poisson equation for the gravitational 
field. Section fVIII deals with the generalization of the Z2 model to a large number of components, partially solving 
the various difficulties pinpointed in the discussion of the most general case. A summary of the results and further 
discussions are presented in the concluding section. 



II. A GENTLE WARM UP: A 2-COMPONENTS SYSTEM 



To clarify the ideas, it is useful to first discuss the simplest instance of a two components systems (similar to the one 
in p^|). in which we have further simplified the dynamics by imposing a 1i symmetry. While this symmetry reduces 
the number of the independent free parameters, it still allows for nontrivial possibilities for the emergent structures. 
It will also help us to elucidate the pattern that we are following in the general case of many components. 

As usual, the system is treated in the second quantization formalism [20(. The quantum field equation reads: 



ft 2 V 2 
2m 



-M + k|*a| 2 + o-|* s | 



(1) 



while the equation for B is obtained by this by the replacement A «-» B. The system obtained in this way is Z2 
symmetric. The masses of the atoms are the same, and as such we have the same chemical potentials, the quartic 
interaction couplings k, a and the U(l) breaking coupling constant A. 

As far as the symmetries are concerned, we see that, of the global U(l) x £7(1) symmetry present if A = 0, 
only a global U(l) survives, corresponding to the invariance over simultaneous phase shifts of the field operators 
$>A,b ~ > exp(ia)^ a,b ■ According to general arguments (Goldstone's theorem), then, we expect at least one gapless 
mode, in the spectrum of the excitations. 

Condensation, as usual, will be treated within a mean field approximation (MFA). Of course, this mean field 
description is just an approximation, and a further analysis is needed in order to test it. In particular, a formalism 
extending the particle number conserving one [U [22J. However, while this is certainly a point to keep in mind, in 
the case of almost homogeneous and weakly time dependent configurations, the MFA does offer a reliable description, 
for ordinary BECs. Clearly, this must be checked with care for the general case we have in mind, but this is a 
development which is beyond the scope of this paper. Within the MFA, the second quantized field operators are split 



1 We will look, in particular, for universality of free fall. 



into a classical part and a residual fluctuation component as 

*x(x) « ?AxI + Xa(x). 



(2) 



Furthermore, to make comparisons with weak field limit of gravitational theories, we will consider the condensate 
wavefunctions to be slight deformations of constant configurations 



V>x(x) = n][ 2 cxp(-i9 x ) (1 + <Ax(x)) , 
where nx,9x are constants and <f>x <C 1. The corresponding Gross-Pitaevski equations read: 

ft 2 V 2 



ihd t ip A (x) 



2m 



fi + K\ip A \ +(t\iI>b\ 



■0a(x) - A^s, 



(3) 



(4) 



and the one with A ^ B. In the limit of homogeneous condensates, these equations become: 

= [-/j + nn A + an B ] n[ 2 exp(-i9 A ) - X^ 2 cxp(-i9 B ) = (5) 

and the same for the other one, obtained by the replacement A ^ B. Let us look at the Z 2 symmetric solution: 

n A = n B = n, A = 9 B = 0, (6) 

where we have used the global U(l) symmetry to fix the phases to be equal to zero. The equation for n, then, gives 
as a result: 



u + X 



k + a 



(7) 



The spectrum of quasi-particles around this configuration can be easily derived from the structure of the equation 
of the fields Xx, in the limit in which we neglect backreaction effects on the condensate: 



h 2 V 2 
2m 



— /j, + 2nn + cm 



\ A + nnx A - \XB + crnx B + cmx B - 



(8) 



By defining the isovector v = (^a A (k), a B (k), a A (—k), a B (—k)j T the equation of motion for excitations can be put 
the matrix form: 



where 



and 



fi 2 V 2 

ihd t v = ( — ^- J + 1U) I ''■ 



10 
10 
0-10 
-1 



(2k + <r)n — [i —A + an nn an 

—A + an (2k + a)n — u an nn 

—Kn —an — [(2k + a)n — a] A — an 

—an —Kn A — an — [(2k + a)n — fj] ) 



(9) 



(10) 



(11) 



With a Bogoliubov transformation from the ax to new operators bx (such that the new operators satisfy the bosonic 
commutation relations), the equation of motion is diagonalized to the form 



ihdtv = 



/wi(k) \ 

uj 2 (k) 

-wi(k) 

V -w 2 (k)/ 



v. 



(12) 
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The eigenfrequencies can be easily computed with this method. They are: 



wx(k) 



K 2 k 2 

2m 



+ 20+ A) 



H 2 k 2 
2m 



1/2 



(from which we see that stability requires fi + A > 0, condition that is equivalent to the positivity of n) and 



w 2 (k) 



h 2 k 2 
2m 



k(3A + m) + {X- v)<7 h 2 k 2 



k + a 



2m 



4A 



k(2A + mu) — /icr 
(k + ct) 



-,1/2 



(13) 



(14) 



We see that the second mode has a nonzero gap, and therefore will represent the massive boson we are interested 
in. The massless mode is the Goldstone boson associated to the breaking of the global U(l) symmetry changing 
both the phases of the field operators "J. Notice also that the gap is proportional to A. In fact, this mode is the 
pseudo- Goldstone boson associated to the soft breaking of the global £7(1) symmetry, which is broken precisely by 
the coupling term (parametrized by A) between the two species. 

If we want to recover the special relativistic spectrum, at least at low energies, we need the coefficients of the 
quadratic terms in k in the two expressions above to coincide: 

(k + ct)(M + A) = k(3\ + (j,) + (X-fi)a, (15) 

leading to 



A = — fi. 

K 



(16) 



Notice that a — leads to A — which is another condition in which we have monometricity, but is also the case in 
which the two condensates are perfectly decoupled. 
In the case where (fTo| is satisfied, the two spectra are: 

1/2 



wi(k) 



h 2 k 2 

2m 



2(l + x)(i 



h 2 k 2 

2m 



w 2 (k) 



h 2 k 2 Y , h 2 k 2 2 

+ 2(l + i)ju— \-4jj, 2 x 



2m 



where we have introduced the dimensionless ratio: 



The healing scale is computed in the usual way |2C 



2m 



K 

It turns out to be: 
1 H 2 



1/2 



J heal 



(17) 



(18) 



(19) 



(20) 



1 + x 4m/i 

Even more important than the spectrum, at least in this case, is the exact form of the Bogoliubov transformation. 
In principle, the calculation can be carried on with the full k— dependence of the various quantities. However, the 
detailed discussion of these transformations lays beyond the limits of this paper. The only thing that we need to 
remember about them is that they are nonlocal field transformations defining the quasi-particle fields (associated to 
the physically propagating collective modes) in terms of the fluctuation operators xx ■ 

In order to understand the coupling between the phonons and the "gravitational field" we need to write the field 
equation for xx when the background condensate is slightly inhomogeneous: 

h 2 W 2 1 

XA + k^aXa + ^Xb + (ti/ja^bXb + (ti/jaiPbXb, (21) 



ihdtXA = 



2m 



fi + 2n\tp A \ 2 + a\tpB\" 



(22) 



where now 

ipx = n^ 2 (l + Ux + %)■ 
We see that the terms coupling the backgrounds ux , vx to the quasi-particles are: 

2n{nu A + (jub)xa + 2nn{u A + iv A )\ A + o-n((u A + u B ) + i(v A + v B ))x B + an [{ u A + Ub) + i{v A - v b )]xb- (23) 

From this expression one easily understand that the effect of the inhomogeneities, from the point of view of quasi- 
particles, is not only a scattering off a potential, but it can also be a transition to a different color. Instead of looking 
for more details, which will be carefully examined later, we will pass to the analysis of the potentials. 
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A. The potentials 

As we have seen before, the quasi-particles will be coupled to the inhomogeneities of the condensates in a variety 
of ways. It is important to distinguish among them, the one that is going to play the role of the "Newtonian 
potential" . Therefore we need to understand what happens to the Gross-Pitaevksi equation when we insert explicitly 
nonhomogeneities. It is immediate to see that the equation will have the following form: 







ft 2 V 2 
2 m 



fi + (2k + <j)n 



4>A + {en — X)(j>B + Kti(j)A + cmcpB + source, (24) 



where source denotes various expectation values generalizing the anomalous mass and the anomalous densities of the 
single BEC. Taking the real and imaginary parts, we can put the equation in the form: 



V 



n 2 v 2 

2m 



/(2 


+ x) 


X 







( UA 




X 


(2 + x) 





=•)) 


U B 


V 








X 




VA 








—x 




\ V B 



source. (25) 



This matrix is block diagonal. While the modes arising from the it-fields are short range, among the ones arising 
from the v— fields there is one that is long range, as a direct calculation can show. In particular, the mode: 

Vl on g = v A +v B (26) 

will obey an equation of motion of the form: 

h 2 V 2 

— Viong = source, (27) 

2m 

which has to be compared with the Poisson equation for the gravitational field: 

V 2 $at = inGpmatter- (28) 

The correct mapping between Vi ong and $jv is just a matter of keeping track of all the dimensionfull coefficients 
appearing in the equations (|27p and (|2"3")1 (see [ljj for the single component case) . In conclusion, in this simple model, 
there is room for emerging Lorentz invariance at low energy for all the particles, while there is simultaneously the 
possibility to coupling the massive mode to a long range nonrelativistic potential. 

We have used the most symmetric Lagrangian that allows for at least one massive mode. One could have obtained 
monomctricity without tuning of parameter by enhancing the symmetry from Z2 to SU(2) by setting a = k and 
A = 0. The price to pay is that in this model, as one can see, there are no gapped modes, and therefore this extended 
symmetry excludes the possibility of discussing the nonrelativistic limit of a massive relativistic particle. 

Notice also that the presence of a symmetry in the fundamental action is not enough to enforce monomctricity. 
Within all the symmetric Lagrangians compatible with the requirements, one has still has to tune the parameters of 
the model [1 81 ) . This fact will carry over in the general multi-component case that we are going to discuss now. 

III. THE GENERAL MULTI-COMPONENT MODEL 

Having examined in a certain detail the case of a two-component system, we can now give the description of the 
general multi-component case. Every component is associated to a second quantized field operator, 

^« = -^E^( k ) elk ' X ' ( 29 ) 
v k 

where A = 1, ...,N C is labeling the components. For convenience, we will often refer to the index A as a color index, 
even though it is not necessarily associated to a representation of some internal symmetry group. For instance, it 
could be associated to a spin component as in the so-called spinor condensates [23[ or to some more general hyperfine 
level like in two components BECs. 

The Hamiltonian of the system is given as: 



//= / ( ~> y\(x)^*A(x) + V@)Y (30) 
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where vha denotes the mass of the atoms of the A-\h color and V is a potential term controlling the interactions 
among the particles. The general structure of this term is not specified. In the rest of the paper we will try to establish 
a number of conditions which must be imposed onto the potential in order to match the properties of a Newtonian 
theory of gravity. 



A. Condensation 



As usual the condensation is described with the mean field approximation (MFA): 

*a(x)«^(x)I + 6xa(x), (31) 

where tpA represents the condensate wavefunction associated to the color A, \A represents the corresponding fraction 
over the condensate and e <C 1 has been introduced in order to explicitly keep into account the fact that the 
noncondensed fraction is small. It might be thought as a 1/N term, where N is the total number of particles, if one 
imagines a sort of asymptotic expansion in the number of particles. 

In the limit of e — » 0, the condensate and the noncondensate parts can be thought as weakly coupled (in particular 
the "backreaction" of the noncondensate fraction on the condensate can be neglected), and the analysis of some 
physical features of the system are easier to be described. In particular, the equation of motion for the field operator 
becomes an equation of motion for the classical fields ipA- 

b^a h 2 w 2 , dv , ,, 

ih ^r = -^ A + WA W - (32) 

These equations are the generalization of the Gross-Pitaevski equation to the case of the multi-component condensate. 
Among the solutions, we are interested into the constant ones, associated to a flat background acoustic geometry. We 
will denote them by V~A- They solve the equations: 

dV - 

W - 0, (33) 



dip' 

and their complex conjugates. 

The number of these constant solutions depends on the number of the minima of the potential V. Symmetries of 
the potential might also affect the structure of the set of these "vacuum states" . The detailed discussion of the_se 
vacua is immaterial for our purposes (at least at this stage), and it is sufficient to focus on a single configuration ipA 
without digging the full structure of the possible configurations. 

Before passing to the discussion of the quasi-particles, it is useful to give some relations between the various 
derivatives of the potential with respect to the field arguments. In order for the quantum system to be consistent, 
besides the positivity of energy which is ensured by the shape of the potential V, the Hamiltonian must be an 
Hermitian operator. Consequently, the potential must be an Hermitian operator. When translated in terms of a 
complex function of complex variables, the hermiticity constraint can be written as: 

vw,r) = v(r,Tpy. (34) 

Therefore, when taking derivatives with respect to the arguments, we get: 

dV fdVY (35) 

(36) 



d 2 V ( d 2 V \* d 2 V ( d 2 V 



dtpAdtpB \dip\di/>%J dtp A dtp B \dip%dil> B 

These identities will play a role in slightly simplifying the following analysis. 



B. Quasiparticles 

In the expansion in power of e, while the lowest order equations are giving the GP equations for the condensate 
wave- functions, the next order is giving the equation for the propagation of the fields \A- These equations are: 

ih^ = -^XA+X A BXB + Y A BX B , (37) 
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where 

d 2 V d 2 V 
Xab = 7D — k~, — > Yab = tt, — ttj-t (38) 

are constant matrices. They are responsible for mixing between different colors as well as for mixing between particle 
and hole states. The propagating modes, the quasi-particles, must be obtained by a suitable Bogoliubov transformation 
of the original modes into eigenstates of the Hamiltonian. Without entering into the details of the calculation, the 
general relation between the quasi-particle field operators (j) A and the particle held operators can be written in a very 
general way as: 

4>a(x) = J d 3 y( y f A B{x 1 y)x B {y)+gAB{x,y)x { B {y)^ , (39) 

where the kernels Jab , 9ab are calculable directly from the Bogoliubov coefficients usually computed in momentum 
space 2 . The diagonalization procedure leads also to the spectrum of the quasi-particles, 

ih^- = E A (-v 2 )4> A + Ef4> A , (40) 

where the term E^ has been introduced in such a way to explicitly show the energy gap of the mode. In terms of 
the wave number k, the energy of the mode is 

E A (k) = E A (k) + Ef, E A (0) = Ef. (41) 

The presence or the absence of a gap is related to particular symmetries of the Hamiltonian and of the background 
solutions tf)A- For instance, if some global symmetry is spontaneously broken by the choice of 'ijj, then some gapless 
Goldstone bosons will necessarily arise. 

In order to translate the gaps into rest masses, a notion of speed of sound is needed. Therefore, the structure of 
the dispersion relations of the quasi-particles must be investigated in more details. At low momenta, the dispersion 
relation can be expanded in powers of p = h\k\: 

E A (p)^e hA p + e 2 ,AP 2 - (42) 

If the mode is gapped, in order to have a similarity with non-relativistic particles, the coefficients C\^a must vanish: 

ei.A = 0. (43) 

Obviously, these are conditions which are imposed upon tp (and hence indirectly on the potential V). There are as 
many conditions as the number of gapped modes. Furthermore, the second term in the expansion must be read as a 
nonrelativistic kinetic term: 



^ = Wu> (44) 

where AAa, which is defined by this equation, is the effective mass for the quasi-particle A. With these definitions 
one can provide a notion of "speed of sound" 3 by means of the gap energy: 



2 



(0) 



4 0) =^c^^c^ = ^. (45) 

Ma 



Monometricity conditions are then equivalent to the requirements: 

< A = c 2 ,VA (46) 



2 For a detailed discussion of the implications, see |l7l l25(l . 

3 This definition is only formal. The correct procedure should involve the discussion of the group velocity in the relativistic regime. 
However, this is not available for the present discussion. 
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They represent a certain number of constraints on the potential V (and, of course, on the solutions ~tp A ). It is rather 
complicated to disentangle them and write down a simple condition on the potential V. For the purposes of this 
paper it will be enough to leave them in this implicit form. For a detailed discussion of this problem in the case of 
the 2BEC system see 

Here we find the first tuning issue. As in the case of 2BECs, the monometricity of the emergent acoustic structure 
is far from obvious. In particular, there are distinct low-energy geometric phases (monometric, multimetric and 
genuinely Finslerian [III), each of them corresponding to a given region of the parameter space for the microscopic 
parameters. 



IV. BEYOND GROSS-PITAEVSKI 



So far, the discussion involved just the homogeneous solutions, where the backreaction effects due to quasi-particlcs 
are completely discarded. To include them, it is needed to change slightly the formalism and to pass to the Bogoliubov- 
de Gennes one. This is still a MFA method, but now the effect of the noncondensed fraction is systematically 
introduced and taken into account. As for the Gross-Pitaevski equation, there is a straightforward generalization to 
be done in order to include the various colors. A consequence of the mean field approximation is that 

(xa) = 0, (47) 

so the first nonvanishing expectation values are the quadratic ones: 

tTUs(x) = (xa(x)Xs(x)}, n AB (x) = (x a (x) Xb (x)), q AB (x) = (Xa( x )Xb(x)), (48) 

which are generalizing the anomalous mass and density of the single-BEC system. A normal ordering with respect to 
the basis x removes an unphysical zero point energy 4 , as well as terms of the form (xx )• 



A. Condensate 



The general study of the properties of nonhomogeneous condensates is a topic which is beyond the scope of this 
paper. Rather, it will be enough to restrict the analysis to small deviations from homogeneity: 



ipA ~ tpA + u A (x) + iv A (x), 



(49) 



which corresponds, as we have already mentioned, to the weak field regime of the corresponding acoustic geometry. 
The slight discrepancy with the notation used in section [TT] is due to convenience. 

Using the MFA into the equation for the evolution of the operator ^> and taking the expectation value, we obtain an 
equation for the condensate which explicitly takes into account the (lowest order) contribution of the noncondensed 
fraction. If one neglects higher order expectation values, which can be thought to be further suppressed by powers of 
1/yN with respect to the anomalous masses and densities, one finds the following equations: 



d ~h 2 \J 2 

ih—(u A + iv A ) = — (u A + iv A ) + X AB {u B + iv B ) + Y AB (u B - iv B ) + XU, 

at 2m A 

where Y, A = + iT, T A are the source terms encoding the effect of fluctuations around the mean field, 



dX Al 



dip c 



xo.bc 



^p—tjj 



( OXab 



9Yai 



c 



n B c 



dY Al 



1p — 1p 



dip* 



ijj—tp 



(50) 



(51) 



Taking the real and imaginary parts of the Eq. (|50[) , and by considering quasistatic configurations, for which time 
derivatives give sub-leading contributions, the equations for the fields u A and v A are easily obtained 



ft 2 V 2 

2771,4 



U A - il AB U B - ^ AB V B 



yR 



(52) 



4 Notice that this normal ordering is not equivalent to normal ordering with respect to the quasi-particle field operators <j>. 
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fi 2 V 2 



-VA 



Yabvb - Tabub = S^) (53) 



2mA 

where, for convenience, we have introduced the notation 

n AB = Rc(X AB + Y AB ), Sab = lm(X AB - Y AB ), (54) 

Tab = Rc(X AB - Yab), Tab = lm(X AB + Y AB ). (55) 

These matrices control certain mixing terms between the different perturbations of the condensates, which are the 
candidates for representing our emergent interactions. These mixing terms must be reabsorbed by a suitable diago- 
nalization, in such a way that we are left with equations in the form: 

V 2 ua + YT^ A = source - (56) 

A 

However, for this diagonalization to be performed, some conditions must be met. In particular, let us introduce the 
matrix M = diag(mi, ...,mjv c ) and the the two 2N C x 2N C matrices, 

o Mr ( r t ) ' ( 5? ) 



such that the equations can be written as 



h 2 

—V 2 (A.w) + B.w = 0, (58) 



where w = (u\ 1 u^ c , Vi, vn c ) t . 

Then, the condition we are seeking is that A and B must be simultaneously diagonalizable. Since M is already 
diagonal, this poses strong constraints on the matrices S,r, Y,f2, unless the degeneracy of the eigenvalues of M is 
large. If M = ml, then the only remaining condition is that IB is symmetric 5 . This provides three matrix conditions 
on each block: 

n T = 0, T T = T, T T = 5, (59) 

which are ultimately conditions on the second derivatives of the potential V, evaluated on the background solution 
■ipA- If these are not met, then the LHS of the equations cannot be put into diagonal form and there is an unavoidable 
mixing between the various potentials. 

However, it is easy to see that the very definition of the matrices Xab and Yab (and in particular the condition 
that the Hamiltonian is Hermitian) ensures that these conditions are identically satisfied. Indeed, if V is sufficiently 
smooth, 

d 2 V d 2 V 

Xab = » , » . = » , » , = X B a, (60) 
dlpAOIpB dip B dtp A 



which ensures that Xab is a symmetric matrix, while 

8 2 V 

dip A dip* B \dip* A dip B 



d 2 V ( d 2 V , 
Y AB = ^4-^=[^J- =Yba, (61) 



which implies: 

Re(Y AB ) = ^ D - ' AB - '™ ' = Rc(Y BA ), (62) 



Yab + YX B Y AB + Y BA 



t iv \ Yab — YXb Yab — Y B a t ^ ,„ , 

lm(Y AB ) = = = -lm(Y BA ). (63) 

2,1 2i 



We have turned our complex equations into real equations in twice the number of variables. Therefore the matrices we are playing with 
are real. 
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Using these properties it is immediate to prove that VLab = &ba and Tab = T^, from the symmetry of the matrices 
Xab and ~Rb(Yab), while the antisymmetry of the matrix In^Y^) ensures that V ba — ^ab- Therefore, the mixing 
matrix between the various potentials is symmetric, and hence diagonalizable. Consequently, the diagonalizability 
condition 

[A, B] = (64) 

can be seen ultimately as a condition on the degeneracy of the matrix M. 
If the diagonalization can be performed, then the equations take the form: 

V 2 u a + \a^ua^Pa^b+Qa^b, (65) 



2mA 

where p, q are the matrices related to some blocks of the diagonalization matrix, niA are the eigenvalues of the mass 
matrix, and Xa are eigenvalues of the matrix IB. Of course, analogous equations hold for fields va- 

It is the case to spend few words concerning the necessity of this diagonalization. Of course, from the point of view of 
the physics, the potentials it, v have the most immediate relevance, being directly related to fluid dynamical variables 
of the condensate (via the Madelung representation). Nonetheless, u,v represent the modes (of the condensate) that 
can be excited independently (at least at the lowest order in perturbation theory). This is why they are so important. 

The terms Xa are related to length scales through: 

1 2X A m A , . 

U A = -R5— (66) 

In the case of the single BEC, this turned out to be the healing length, which represents the only dynamical scale 
present in the game. This scale is responsible for the extremely short range potential emerging in that example. 
However, having at our disposal more components, and in particular several interaction terms between them, it 
happens that more interesting situations can be realized. While a general rule is impossible to be discussed, it is 
illuminating to show what could happen in a simplified case. Let us assume that for some reason, related to the 
symmetry structure of the system, as in section [ill one ends up with a set of coupled equations like: 

W!0 + (^VlM source, (67) 



then, diagonalizing, one finds: 



U-2 I \ a a I \ U2 



u\\ ( 2a A / u\ 



v U + o o U > =source ' (68) 



where one of the two potentials is short range, while the other is long range, giving rise to 1/r potentials, rather than 
Yukawa ones. Therefore, while the interaction between components complicates the analysis of the physical properties, 
it opens the possibility of creation of long range potentials which are strictly forbidden in a single component BEC. 
However, as one can imagine, this feature is not generic at all. Some requirements are needed concerning the (global) 
symmetries of the system on one hand, and on the symmetries that are preserved by the particular choice of the 
vacuum state on the other. In particular, a situation like the one briefly discussed here will be realized in the system 
considered in section IVHI 



B. Quasiparticles 



To study the properties of quasi-particles, in principle one should go for the Bogoliubov transformation in the case of 
non-homogeneous backgrounds. First of all, rewrite the equation for the excitations in non-homogeneous backgrounds 
as 



d_ 

"St 



X(x) 
X*(x)t 



X(x) 
X(x) f 



■W(x) 



X(x) 

£(x)t 



(69) 



where Q denotes both the kinetic operator and the linear term for the homogeneous background while the IA matrix 
encodes the inhomogeneities alone, which are contained just at the linear order, in the approximation scheme we are 



using. Doing the Fourier transform with respect to the modes e lfcx , we get 



ih 



d 



a{k) 



dt v at(-fc) 



9(k) 



a(k) 
ot(-fc) 



+ / d 3 xe 



-ikx. 



W(x) 



X(x) 
X(x)t 



(70) 
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The last term can be manipulated by introducing the Fourier transforms: 

U(x) = J d 3 he ihx U(h), (71) 

to get, finally 

< ( »% ) = s <*> ( s% ) + / ^ ( It-t', ) ■ < ra » 

This means that the Bogoliubov transformations, in this case, are much more complicated and involve different 
momenta, according to the spectrum of U(x). This is expected, since the simplicity of the Bogoliubov transformation 
for the homogeneous background was deeply related to translational invariance, which is broken by the presence of 
the inhomogeneities. 

Given the intrinsic difficulties in disentangling Eq. (j72[) . we will not try to compute the Bogoliubov transforma- 
tions in this case. Rather, we will try to deduce the Feynman rules for the interaction between phonons and the 
inhomogeneities, using a scattering theory point of view. 

We will use, then, the free states obtained linearizing around the homogeneous background, i.e. the phonons 
described in section MI B[ without including the effect of the corrections to Bogoliubov transformations. We will 
treat these states as asymptotic states, in the sense of scattering theory, imagining that the inhomogeneities are small 
localized ripples in the condensate. This is another limitation of this analysis but, nonetheless, we will be able to grasp 
the essential features of the physics of the system. 

The perturbation in the Hamiltonian for the quasi-particles is encoded in the operator 

H in t = J d 3 x (a AB {x)x A {*)XB(*) + b AB (x)x\(x)x B (x) + h.c. ) , (73) 

where the (matrix) functions a,b are linear combinations of the perturbations u A ,v A . The explicit form of these 
relations are not particularly illuminating, so we will not report them here. 

This operator can be written in terms of the (free) quasi-particle fields by using explicitly the Bogoliubov transfor- 
mations 6 

Xa{x) = J d 3 y (r AB (x 7 y)(p B (y) + s AB (x,y)<f> B {y)J (74) 

and the Hermitian conjugates. It will suffice to look at the first term, since the others can be treated in the same 
way. So: 

His = / d 3 xa AB (x) X A(x-)XB(x) + h.c. = 



^yd 3 z (a A I B (y,z)0 A (y)0B{z) + b A I B (y,z)4> A (y) f 4> B ('z)+ h.c. 



(75) 



where the kernels a, b are obtained by the convolution of the Bogoliubov kernels with the matrices a, b. For example, 
we have: 

<4s(y. z ) = J d 3 x\ i a RS (x)r RA {x,y)r SB (x,z) + (b RS (x)s RA {x,y)s SB {x,z)y^ . (76) 

As said, we will be interested in matrix elements between phonon states of definite momentum. Given the quadratic 
nature of the peculiar "interaction term" that we are considering, the interesting matrix elements (in the Born 
approximation) are those between single quasi-particle states with definite momentum, i.e. 



T(h c ,k D ) = (k D \H mt \h c ). 



(77) 



Actually, these are the inverse of the transformations considered in lHI Bl 
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By straightforward manipulations, we see that 

T(h c ,k D ) = f d 3 yd 3 z(b A B(y,z)(0 gp |&c(k)^(y)^(z)6yh)|0 9f) ) 



+ b 



\ B (y,z)(O q /b c (k)$\(y)$ B (zp D ^)K P )) , (78) 



whence 



T(h c , k D ) = J d 3 yd 3 z (b CD (y, z) e -* z e^ + b f DC (y, z^^e" 12 ) + BS DC S 3 (k - h). (79) 
In this last equation we notice the appearance of a forward term with a coefficient 

B = [ d 3 xS AB b AB (x,x), (80) 



signalling the fact that a wavefunction renormalization must be performed. Notice that this coefficient is crucially 
dependent on the structure of the Bogoliubov coeffients but also on the perturbations ua,Va- 
There are also quasi-particle pairs production/destruction graphs, corresponding to 

T({h c ,k D },0) = (k c ,h D \H int \0), (81) 

for which the other coefficients a would be relevant. However, these graphs are not relevant for the discussion of 
Newtonian gravity, which involves only scatterings off a potential, without particle creation effects (at least not in the 
regimes where one usually applies the Newtonian limit, i.e. well below the threshold for particle production). 

Notice that, in principle, these matrix elements can be opened up and rewritten in terms of suitable integrals of 
u Al vb, which are still entering linearly in these matrix elements. 

The comparison with the standard Newtonian potential requires much more work with respect to the single BEC 
case. We have to compare these results with the ones obtained via 

-^Newtonian = ~ / d 3 X ^ M A ^Newtonian (*) 4> A (x)(f) A (x), (82) 
A 

where Ma are the masses of the quasi-particles. In particular, one has that 

T N (h c ,k D ) = -S CD Mc I d 3 ^Ncwtonian(x)e l ( h - k > x (83) 



T N ({h c , k D }, 0) = T N (0, {h c , k D }) = (84) 

Therefore, in the multi-BEC case, the identification of the Newtonian potential involves the comparison of all the 
expectation values given above with the ones given in the Newtonian theory, leading to a much less trivial set of 
conditions. In particular, we want that 

T(C, D) s» T N (C, D) ocM, (85) 

where M. is the mass matrix for the quasi-particles. As one can easily see, this amounts to ask N c (N c + l)/2 conditions 
on the matrix elements (f78| (for each pair of momenta!), which are ultimately related to the shape of the potential 
V describing the condensate and to the homogeneous solution playing the role of our background. 

Since T(C, D) is a Hermitian matrix, and we can extract its projection along the mass matrix Ai in such a way to 
identify uniquely the candidate for the analogue of the Newtonian potential: 

^^<b ,VW(h-k)x _ tr(T sub {k A ,h B ).M) 

tr(Al) 

where T su b is the transition amplitude without the forward term. 

This way of manipulating the various quantities shows that the (would be) Newtonian potential is a complicated 
linear combination of the various potentials u, v. Furthermore, given that the the relation (|86[) is actually momentum 
dependent, the relation between the Newtonian potential and the potentials u, v is nonlocal, when translated in real 
space, with the following general shape 



^Newtonian 

CD 



This discussion shows how complicated is the identification of the analogue of a single Newtonian potential in the 
multi-BEC case, among the various possible candidates (after all, we should expect from the very beginning a sort of 
multi-metric structure, whence a multi-potential theory). 
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V. SOURCE TERM 



So far the source term has been neglected. However, it is interesting to spend some time to study it, because it 
will lead to additional properties of the emergent theory. As we have already said, the source terms appearing in the 
Bogoliubov-de Gennes equations have the following structure, at the lowest order: 

„ ( dX AB \ ( dX AB dY AB \ J<> y ab\ .... 

where all the expectation values are written in terms of the helds x- 

As in the case of the single BEC, these terms contain vacuum contributions related to the depletion factor pol |24|. 
which however has a more rich structure, in this case. Inserting the Bogoliubov transformations, we get that the 
expectation values in the Fock vacuum for quasi-particles are: 

m vacuu m{x) ^ JePyeP Z r AC (x,y)s% D {x,z){$c{v)&D{*)h (89) 

whence, using the commutation properties between the fields <f>, one gets: 

m vacuum {x)= f d 3 y J2r AC (x,y)s BC (x,y). (90) 
J c 

Similarly: 

rffi mm {x)= [ d 3 yJ2sAc(x,y>Bc(x,y), q%g"> m (*) f d 3 y £ r* AC { Xl y)s BC {x, y). (91) 
^ c c 

These terms are generating shifts in all the potentials due to vacuum effects (the inequivalence between the Fock 
vacuum for the atoms and the one for the quasi-particles). In some sense these terms migh be interesting because 
they could be related to small cosmological constant, as suggested in [13] ■ However, here we will neglect the effect of 
quantum vacuum fluctuations and consider the effect of quasi-particles. 

Coming back to the source term, it is obvious that the interesting parts involve states with some quasi-particles. 
The expectation values to be computed are of the form: 

m = (n\xAXB\n) , (92) 

where |n) is a stenography for an n— quasi-particle state. Similar expressions hold for n, q. To be concrete, let us 
study, as in [13], the case of a state like: 

\()=J2 [d 3 x( A (x)4> A (x)\n). (93) 

A J 

This state corresponds to some quasi-particles which have a certain number density specified by the "wavefunctions" 
Ca(x) (we are using the normalization ((\() = 1). Denoting with n A (x) the number density of quasi-particles of color 
A, it is easy to verify that 

n yl (x) = (C|^(x)0 A (x)|C) = |a(x)| 2 . (94) 

In order to reproduce the correct Newtonian behavior, we need that the source terms become function(al)s of these 
number densities. However, this is not the case. Consider, for instance, the term: 

m AB [C] - (CI£4(x)xs(x)|C). (95) 

Using the definition (|93p . the Bogoliubov transformations and the commutation relations between the quasi-particle 
fields <j>,(fy, one obtains: 

m AB [Q] = m A T um +J2[ d 3 yd 3 zA ABC D(x;y,z)Cc(yKD(z), (96) 

CD J 
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where the kernel A, related to the Bogoliubov coefficients, is: 

AAsci)(x;y,z) = ^ D (x,z)s BC (x,y) + sac(x, y)r SD (x, z). (97) 

Similar results are obtained upon examination of the other terms obtained from n, q. The most important lesson 
that must be taken from these computations is that the source term is a complicated nonlocal functional of the 
quasi-particle distribution, not in terms of the number density of quasi-particles (as one naively would expect), but 
rather in terms of their quantum wavefunction (|93p . 

While this is of course a crucial difficulty in obtaining some meaningful source term related directly to the number 
density, this feature should not come as unexpected. In fact, quasi-particles are collective motions: this represents 
an interesting example, albeit very simple, of mismatch between microlocality and macrolocality pH ]. Without 
any additional ingredient able to hide it, the equations of motion for the effective theory of perturbations around 
homogeneous background will be displaying explicitly nonlocal effects. This feature was already present in the case 
of the single BEC. In the case of several components, it is easy to see that the effect is worsened by the appearance of 
off-diagonal terms C/t(y)C-B( z )j with A ^ B, which do not have any sort of interpretation in terms of single components. 



VI. EMERGENT GRAVITATIONAL DYNAMICS? 



Up to now we have already listed a number of difficulties that one has to face in order to build successfully a 
model of emergent Newtonian gravity starting from a multi-component BEC. Yet, we do not have addressed a most 
important point: the equation of motion for the "gravitational" field. 

In section IIV1 examining the properties of the matrix T encoding the scattering amplitudes of the phonons off the 
inhomogeneities, we have seen where the "gravitational" potential is hidden, i.e. in the part proportional to the mass 
matrix of the quasi-particles. From our perspective, to recover something resembling Newtonian gravity, we cannot 
only require that all the remaining matrix elements of T must be small, but also that this "Newtonian potential" 
must obey a Poisson-like equation, 



1 

L2 



^Newtonian 



with L a constant determining the range of the interaction. Since ^Newtonian is a nonlocal linear combination of the 
potentials u, v, as we have previously mentioned, the equation of motion for ^Newtonian will be obtained by considering 
the equations of motion for these potentials. We have established that: 



A, 



^NewtonianOc) = / d 3 y ^ (a A (x, y)u A (y) + Pa (x, y)v A (y)) , (99) 



A=l 



where the precise structure of the kernels a,j3 can be computed from |87|) . Taking the Laplacian of this expression 7 , 
we obtain that: 

- Ac 

V 2 $Ncwtonia„(x) = d 3 yJ2 (V^ A (x, y)«U(y) + V2/3 A (x, y)v A (y)) , (100) 

J A=l 

where we have explicitely introduced the notation V 2 to remember that, when acting on the kernels a, (3, the operator 
should act on the x variable. This result shows that in order to get to some form of Poisson-like equation we have to 
unveil the structure of the kernels a,j3, which, in turn, are determined by equation (|87j) . 
By comparing (jlOOp with the desired result, 

V 2 <f>Ncwtonian(x) = ^Newtonian + C A 12 A e + d A Yj A m , (101) 

we see that the kernels must not be generic. They must obey some integral identity of the form 

J d 3 yVla A (x,y)f = t a J d 3 ya A (x,y)f + c A J d 3 ya A (x,y)V 2 y f. (102) 



7 We assume also that the kernels are smooth enough to be allowed to derive under the integral. 
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Albeit already very stringent, this condition is not enough to ensure that the equation will have a structure like 
the one we want. To see this, let us assume, for simplicity, that the kernels are proportional to Dirac distributions, 
i.e. that the Newtonian potential is a simple linear combination of the potentials u, v. 

^Newtonian 

(x) = (a A UA(x) + PavaW) ■ (103) 

A=l 

In this case, the action of the Laplacian is manifest: 

V 2 $Ncwto„ia„(:c) = J2 ("aV 2 u a (x) + ^aV 2 w a (x)) , (104) 

A=l 

whence, using equations (|56|) (and analogous for the v), we get: 



V 2 $Nowtonian(a:) = ^ a Aj2 *U(x) + Pajq, ^( x ) + source. (105) 

A=1 V L A,(u) h A,{v) J 

However, unless the La,u,La,v are all the same, we are not able to reconstruct the needed ^Newtonian to obtain a 
Yukawa-mass- like term. This implies that, even if we discard the corrections due to the intrinsic nonlocality of the 
Newtonian potential, we have to tune the constants determining the La for the ua,vb, in order to make them to 
have the same value. 

Besides the important difficulties concerning the left hand side of the Poisson equation, another important issue 
involves the emergence of a unique Gn, or, equivalently of some sort of universal law of gravitation. To understand 
whether this can happen, the equations for the potentials must be reconsidered with care, and in particular the sign 
of the coefficients that will lead to the source terms for the Poisson equation. 

As it is easy to understand, by doing all the diagonalizations/Bogoliubov transformations etc. we cannot control 
and ensure that all the coefficients in front of the pseudo-densities of quasi-particles have the same value (enforcing 
a universal law of gravitation). In fact, as we have already noticed, we do not have any control on the sign of these 
coefficients (i.e. on the sign of for each component). In other words, we could have some form of antigravitation, 
at least in principle. Furthermore, as we have noticed in the previous section, we can have in the source term some 
sort of cross-terms which depend simultaneously on the distributions of quasi-particles of different colors. 



VII. THE S Nc CASE 



The results of the previous section seem to be pretty discouraging, confirming the naive intuition that without 
Lorentz invariance and without any form of equivalence principle imposed a priori on the microscopic theory, the 
effective theory of phonons will not display these properties at all. 

In this section we will discuss an example in which the presence of an internal discrete symmetry ameliorates 
the situation, leading quite naturally to a form of emergent equivalence principle. Lorentz invariance, as usual, is 
not completely under control, and a tuning of the microscopic parameters has still to be performed. Nonetheless, in 
absence of this tuning, the excitations will be divided into multiplets which are homogeneous with respect to spacetime 
symmetries and to the coupling to what will turn out to be the only possible long range potential in the model. 

The system we are going to discuss is the immediate generalization of the Z2 model that we have considered in 
section |TI] to the case of N c components. The Hamiltonian for the system will be: 

ihd^A = - m) *a + «i*a| 2 *a + [a J2 i*si 2 - A *B- ( 106 ) 

V m / y b^A J B^A 

where A, B = l,...,N c are labels for the bosonic components and /x, k,<t, A the only four free parameters. This 
Hamiltonian is invariant under the full permutation group of the components, 5jv c , and under global U(l) symmetry 
^a — > e ie ^>A,yA. This Hamiltonian allows the largest symmetry group leading to quasi-particles having dispersion 
relations with relativistic regimes. Had we imposed a U(N C ) symmetry, this would not have been the case. 
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The Gross-Pitaevski equations read: 

ihd t ^ A = f-^p - A tpA + k\iPa\ 2 iPa + I o ]T \4b\ 2 J - A ]T 4b, (107) 
and the Sn c -invariant constant solution is 

V'a=" 1/2 , VA, (108) 

where n is the number density of the condensates and the phases are set to zero, by using the global U{\) symmetry. 
Plugging this ansatz in the GP equation we get 



n^ + ^-y. (109) 
K + (N c — l)cr V ' 



A. Quasiparticles 



The quasi-particle spectrum has to be obtained applying the method of Bogoliubov transformations to the particular 
case. Splitting the operator ^a using the mean field approximation, 



i> A = n 1/2 +XA 



(no) 



The linearized field equation for x reads 

h 2 v 2 



2m 



A* j Xa + Zkhxa + Knx A + {N c - l)an XA + an \ ^ ( X b + Xb)\ - A X! XB ' ( m ) 



.B^A 



B^A 



Going in momentum space, and defining the isovector _Fjv c (fc) = ^ai, ■ • • , ajv cJ o,\," • , ajy J T j then Equation (|111|) 
can be written in matrix form: 



/ 



ihd t F Nc (k) = 



A : B 



\ 



F. 



(112) 



V -B : -A ) 



where we have defined the matrices 



A(k) = a(k)l Na +bl Na ; B = cl Nc +dl Nc , 



(H3) 



with 1n c being the N c x N c identity matrix, ljv e the N c x N c matrix having all the entries equal to one, and the four 
numerical coefficients given by: 



a(k) 



h 2 k 2 
2m 



/i + 2nn + (N c — 2) an + A, b = an — A, c = (k — a)n, d = an. 



(114) 



The matrix 1 is real, symmetric and hence can be diagonalized by means of an orthogonal transformation R. It is 
easy to see that its diagonal form is 



( N c ... 0\ 




R 1R= J 



V 



(115) 



/ 



where the eigenvector associated to the nontrivial eigenvalue is (1, 1, 1, 1), when written in term of the original 
basis. 
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As a consequence, we can easily put in diagonal form (jlll|) . First of all, one has to diagonalize the single blocks of 
the dynamical matrix by means of the transformation R: 



[o : R T J 



A 



B 



\-B : -A) 



( R °) 

V o ; r J 



( 



\-D B 



D, 



\ 



Da j 



where 



D A = a(k)l + bJ, D B = cl + dJ, 



(116) 



(117) 



are diagonal matrices. Finally, by permutations of rows and columns, the dynamical matrix can be put into a block- 
diagonal form, where the blocks on the diagonal are two by two matrices having the same structure of the dynamical 
matrix of the single component BEC, namely: 



ai(k) Pi 
-ft -Ot(fc) 



(118) 



and hence leading to the same steps to complete the Bogoliubov transformations. We leave the details aside, and we 
focus on the emerging spectrum. The eigenvalues of (I118P are easyly obtained. The spectrum of the mode is given 
by: 



(fc) = a 2 - ft. 



(119) 



In view of the structure of the matrix J we can conclude that the spectrum of quasi-particles will consist of two 
families. 



• A gapless singlet, when a(k) = a{k) + N c b, f3 — c + N c d, with dispersion relation given by 

h 4 k 4 , „,K 2 k 2 



Am 2 



(w + N c b)- 



ni 



• A multiplet of gapped modes, with dispersion relation given by: 

, h 4 k 4 h 2 k 2 



w- 



Am 2 to 

all with the same rest mass, determined by ujq (and by w, controlling the speed of sound). 
Here, we have introduced 

w = -[i + 2k77 + (N c - 2)0-77, + A 



(120) 



(121) 



(122) 



to simplify slightly the expressions and highlight the relevant contributions. Clearly, in order to have the emergence 
of low energy Lorentz invariance, we have to impose a constraint on the coefficients of the terms quadratic in k in the 
two dispersion relations. In particular, the monometricity condition reads 



6 = an = A. 

This condition leads straightforwardly to a number of simplification. First of all, using f|109[) : 

Xn 

a = — = xk, 



(123) 



(124) 



where we have introduced the dimensionless ratio x = X/ u = a/ k. Then, plugging this result back in (|109[) we have: 

77=^. (125) 

K 

Notice that all these results match with what we have already seen in the Z2 case, in section HU 

It is worth stressing again that monometricity is not imposed automatically by the underlying symmetries of the 
system, but it must inserted by hand with a suitable tuning of the parameters. However, the modes are organized in 
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multiplets (with respect to the residual symmetry Sjv c -i) within which the low energy emergent metric structure is 
the same. Lorentz violating effects, in this kind of scenarios, will be present only when confronting different multiplets. 
This is perfectly analogous to the way in which Lorentz violating terms percolate in a field theory with spontaneous 
symmetry breaking [27[, and confirms the expectation that internal (custodial) symmetries might ameliorate the 
naturalness problems of Lorentz violating theories [HI, [2j| . 

Before going on, we need to introduce the healing length. The healing length should be computed by examining 
the dynamical properties of the condensate (20| . Here we define it by means of the dispersion relation of the particles, 
as the scale at which we have the crossover between the low energy Lorentz invariant theory to the high enery Galilei 
invariant theory, and it is the UV scale at which the low energy hydrodynamic description of the condensate breaks 
down 8 . It turns out to be: 

2 _ h 2 _ / l \ h 2 

h ~ 4m( M + (N c - i)A) ~ ^l + (JV c -l)a:J 4^? (126) 
This scale will set the characteristic scale with which we will make all the comparisons, in the next section. 

B. The potentials 

Let us now pass to the examination of the equation of motion for the condensates, when we turn on the backreaction 
effects. In the quasi-static, weak field limit, the Bogoliubov-de Gennes equation for the condensate turns out to be: 



2V72 



2m 



fx) (u A + Wa) + m(3wa + Wa) + (N c - l)\(u A + Wa) + A ^ (ub - Wb) = S A , (127) 



where E A are given by expectation values of the fluctuation operators, as seen in the previous section. We can put 
this equation in a diagonal form, by taking real and imaginary parts. Defining U = (iti, Ujv e ) T , V — (i>i, hjv c ) t 
and £ = (Si, Sjv e ) T , the equation can be written in matrix form: 

.^-I Nc +aI Nc + bl Nc ^U = Re{T,), (128) 
?L 2 V 2 „ 



2m 



<& No +dl NB )V = 1m(E), (129) 



where 

a = 2/i + (N c - 2)A, b = A, c = N C X, d = -A. (130) 

The range of the effective potentials can be understood immediately by diagonalizing these matrices, so that the 
equations for the effective potentials decouple and get the form: 

-V 2 + j Hi = source (131) 

Using the same technology introduced in the previous subsection, we can easily determine the way in which the 
condensate reacts to the presence of quasi-particles, by looking at the properties of the potentials u and v that 
diagonalize (fT28|) and (fT29]l . 



a. There is a potential u\ coming from the U which has range given by: 

1 _ 2m 

Zf = ~¥ 

Trivial manipulations will show that L\ = L] x 



— = —( a + N cb)- (132) 



The contribution of the quantum potential to the hydrodynamic equations obtained from the GP equation with the Madelung repre- 
sentation becomes relevant when the inhomogencities in the condensate wavefunction reach the healing scale. 
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b. There are (N c — 1) potentials u with a different length scale. The relation with the healing length is easy to 
obtain: 



Ll_ 1 + ^x 



L\ l+(N c -l)x' 



In the limit 9 N c x — > oo, 



Ll 2' 

which signals that these potentials will be associated to short range phenomena as well, 
c. There are {N c — 1) potentials v with range given by: 

1 1 N r x 1 



(133) 



(134) 



L\ 2 1 + (N c - l)x L\ ' 
still associated to short range deformations, as one can see by inspection of the limit N c x — > oo. 



(135) 



d. Most importantly, there is a long range potential, v\, due to the fact that one eigenvalue of the constant matrix 
in (jl29[) is given by 



oc c - N c d = N c \ - N C X = 0. 

Li 



(136) 



Therefore, in this class of models, there will be always a long range mode, coming from the combination of the 
perturbations for the velocity potentials of the condensates 



1 Nc 
N ^ 



v B . 



(137) 



B=l 



C. Coupling to the long range potential 



We want now to understand what is the coupling between the phonons and the potentials that we have discovered 
so far. Instead of using the method of scattering matrix as in the general case, here we focus on the coupling with 
the unique long range mode that we have, setting all the other modes to zero. Trivial manipulations show that this 
corresponds to the particular configuration in which: 



va(x) = $(x), 



VA = 1,...,N C , 



(138) 



with the ua set to zero. We have now to reconsider the derivation of (jllip with our hypothesis of Lorentz invariance 
of the low energy branch of the spectrum and of small inhomogeneities parametrized by (|138p . One gets: 



ihxj 



fi 2 V 2 
2m 



■H ) xa + 2fiXA + ^x\ + 2«M$X^ + 2U$ 

B=£A 



Xb 



(139) 



The coupling of the modes with the field $ can be understood very easily, by making use of all the algebraic tricks 
that we have seen by now. We can rewrite the last equation in matrix form, as before. The dynamical matrix will 
have the form 



' A : B + iC 
\-B-iC : -A ) 



(140) 



This limit has a simple interpretation. Given that the mass of the massive modes is a monotonically increasing function of A, this limit 
corresponds models with a large number of heavy fields. 
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where A, B are defined as in the homogeneous case, and 

C = 2*(x)(0*-A)I + Al). (141) 

Therefore, before the last Bogoliubov transformation, the dynamical matrix can be written in a block diagonal form 
in which each block has the following form: 



(142) 



where the value of qj depends on which eigenvector (and hence which quasi-particle) we are considering. It is clear, 
then, that the modes will be coupled to $ in a non-universal way. The gapless singlet will have a value of the charge 
q — 2(/i + (N c — 1)A which will be different from the value of q = 2/j, for the massive quasi-particles. 

If we want to interpret $ as a form of Newtonian gravitational potential, this is an instance in which even after the 
imposition of Lorentz invariance, the equivalence principle is still violated. This is of course an expected result, even 
considering that one of the quasi-particle species is gapless. However, if we limit ourselves to just the gapped modes, 
which are in fact the ones relevant for the discussion of a Newtonian limit of gravity, within this sector all the modes 
are coupled in the same way to the potential, and hence they will fall in the same way under its influence. 

In fact, we can push the analysis a bit forward and consider also how the quasi-particles are source for these 
potentials. First of all, let us consider the two classes of vevs that we have introduced in section V (the third is not 
there, in this case). 

m A B = (xaXb), n A B = (x a Xb)- (143) 

By means of the rotation matrix R, as we have seen, we can define new fields £ 

£i(x) = (i? T )f xb, (144) 

for which the equation of motion is block diagonal, and for which the Bogoliubov transformation can be done sepa- 
rately: 

|a(x) = J d 3 k (a A (k)b A (k)e^ + ^(fc)SJ 1 (k)e- ik - x ) , (145) 

where b, are quasi-particles creation and annihilation operators. As a consequence of the separation in blocks, we 
have: 

m AB = J2 (R T )a (R T )B(ici D ) = J2 ( rT )a(R T )b f d 3 ka D (k)p D (k)(n« D p -(k))6 CD (146) 

C,D C,D J 

where (n q D p '(k)) is the vev of the number of quasi-particles of (effective) color D with momentum k. Notice that in 
this expression we have the same nonlocal behavior that emerged in the analysis presented in [ItJ , due to the presence 
of the scale dependent Bogoliubov coefficients a(k),/3(k). We do not need to compute the form of xi A b, given that 
we are interested in $, which is a linear combination of Uj, whose equation involves the imaginary part of the source 
term. It is immediate to get: 

,VJV2 -$ = 5> T )f3(lM) (147) 



2m 



Now, in this particular case, we know that: 
whence we obtain: 



L A 



(R T )t = VA (149) 



which implicitly tell us that the quasi-particles will be a source for the field < I > (x) in a universal way. When we 
introduce all the required dimensionfull quantities to establish the relation between the gravitational field and the 
field $(x), we will obtain that the analogue of the Newton's constant Gjv is particle independent, at least in the 
massive sector. 
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VIII. CONCLUSIONS AND OUTLOOK 

The analysis presented in this paper is a (partial) assessment of the idea that gravity is a hydrodynamic description of 
some pre- geometric degrees of freedom, at least for what can be inferred by considering a BEC analogue model 10 . This 
rather long and detailed discussion shows that getting a gravitational theory as simple as Newtonian gravity within a 
BEC analogue model, without a strong guiding principle like gauge symmetry or a suitable custodial symmetry, is a 
very long shot. We have shown that, exactly as in other instances of emergent gravity [3lj . an underlying symmetry 
principle is needed in order to tame a number of features that would be in plain contrast with current observations. 

Summarizing, the emergence of a Newtonian gravitational dynamics (and low energy Lorentz invariance) out of a 
multi-component BEC-like system requires a number of conditions and display a number of nonconventional features: 

• A number of "monometricity" conditions (|43l) . (|46|) are needed to translate the gap in terms of mass, and also 
to guarantee a form of low energy Lorentz invariance. 

• Some diagonalizability conditions (|64[) are required to avoid mixing between potentials. 

• The "non gravitational interactions", i.e. the part of the matrix T(A,B) (|75|) which is not proportional to the 
mass matrix of quasi-particles, must be small (in other words we need to screen the color of the components). 

• To obtain a sensible Poisson equation we need a precise tuning of the parameters to control the ranges of the 
various potentials involved (see discussion after Eq. dlOSp '). 

• In order to get a single Newton's constant we need additional tuning. 

• Nonlocality is heavily affecting the gravitational sector: given the mismatch between particle's and quasi- 
particle's notion of locality, the effective equations of motion are displaying nonlocalities at all the possible 
levels. 

These are, in extreme summary, the results of sections IIIIHVI1 Besides them, we have also shown that under certain 
conditions we can have some interesting situations in which things seem to work a bit better, by considering a model 
in which an additional internal (discrete) symmetry is imposed on the system. This leads to a series of nice features: 

• (partial) taming of Lorentz violation: particles are organized in multiplets, when their dispersion relations are 
considered; 

• a long range "gravitational potential" appears; 

• the quasi-particles are coupled to this gravitational potential in a pattern according to which the equivalence 
principle (in the sense of universality of free fall) is implemented within multiplets. 

With respect to these results, the case of the single component BEC seems peculiar, an oversimplified case just 
as the single field normal mode analysis 32] is too simple with respect to the multi-field case to get an emergent 
spacetime [16]: many of the phenomena proper of multicomponent system cannot be displayed because there are 
too few components. However, the multi-component case has also the potentiality of generating effective long range 
interactions, and hence still the analysis is worth the pain. 

The case of multicomponent systems shows that in order to produce a phenomenologically viable emergent gravity 
scenario, insisting in abandoning the basic notions of geometry Lorentzian geometry, together with the principles that 
govern the coupling to matter, one has to assume a lot of additional conditions. Stated in a more explicit way, without 
having some alternative form of equivalence principle from the very beginning, we should not expect it to emerge at 
the end. In fact, as it has been shown in [31| . it might be that there is an accidental emergent symmetry that takes 
care of the equivalence principle. However, in absence of a detailed full-fledged mechanism of emergent gravity from 
a pre-geometric theory it is impossible to address in a complete way the issue. 11 Of course, this difficulty with the 
equivalence principle is nothing else than the dynamical counterpart of the necessity of a custodial symmetry taking 
care of Lorentz violating terms at low energy in Lorentz violating scenarios (for a recent discussion of this point and 
references, see [H|, sec. Ill C). 



A vaguely similar situation was considered in |30H . albeit the setting is pretty different, both conceptually and technically. 

For a discussion on the limits on emergent gravity models (and additional references) we refer to [33f | and to [34l |. where the role of 

symmetries is discussed with care. 
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From the discussion of the Sn c model in section IVII1 it appears that even the imposition of an additional internal 
symmetry is not able to keep under control completely the system: to enforce Lorentz symmetry at low energy we 
need to perform additional tunings. Furthermore, given that the gravitational field, in the model, is completely 
disconnected from any form of gauge principle based on the Lorentz group, even imposing Lorentz invariance is not 
enough to enforce an equivalence principle. From the point of view of real world, this means that if the standard 
model and general relativity are the low energy limits of a Lorentz violating theory, the underlying theory is not 
generic at all, but must have very specific symmetry properties 12 . 

It must be said that this point of view is only partial, since it seems to be associated to the quantum/emergent 
gravity scenarios most commonly encountered in the analogue models literature, which are based, in one form or the 
other, on some condensation mechanism in a pre-existing space. It might that different classes of phenomena are 
needed, i.e. that one should look for a completely different mechanism for the emergence of spacetime and matter 
fields, based on totally different concepts. In this direction, it is important to mention the cases of matrix models for 
2D quantum gravity [3{| and for higher dimensional spacetimes as well as the Group Field Theories (GFT) [42[ 
recently proposed to implement in a background independent way the functional integral approach to gravity. 

In particular, as it is suggested by matrix models, matter degrees of freedom and geometry could be encoded in the 
same objects, which are pre-geometric and are not obviously associated to some fields living on some manifold. The 
geometric and matter degrees of freedom will appear as we know them only when the statistical system defined by 
matrices is forced to approach criticality [39( . A similar situation occurrs in the case of quantum graphity models [41| . 
The way in which matter fields and geometry (and their reciprocal relationship) are emerging from a pre-geometric 
theory has been investigated also within GFT [43jj. In this case, however, the correspondence with the continuum 
semiclassical limit represented by general relativity is much less understood, even though some progresses concerning 



the appearance of a hydrodynamic limit have been done 44J. In all these scenarios, the very nature of gravity and 



matter fields as excitations of the same underlying structure might automatically guarantee the emergence of low 
energy local Lorentz invariance and of a gravitational theory satisfying some form of equivalence principle. Of course, 
more work is needed to understand whether this is the case. 

Lacking a solid derivation of general relativity from consistent pre-geometric models, we can still gather insights 
from analogue models. In particular, the work presented in this paper is not only elucidating the conditions under 
which we can have, from an emergent systems, phenomenology which is compatible with observations, but also what 
are the potential candidates to detect that spacetime geometry is an emergent phenomenon, in this case Lorentz 
violation, nonlocal phenomena in the various sectors of the theory and possible violation of equivalence principle 13 . 

In this paper we have considered the very simple case of nonrelativistic scalar bosons. However, it is natural to 
ask what happens in a system with several fermionic fields, and under which condition we can have the emergence of 
gauge interactions, besides the gravitational potential. A model like this would be much closer to the kind of physics 
we are familiar with from high energy physics, and allow us to discuss even additional features like the appearance of 
chirality. Of course, this requires us to change completely the class of models (and look for Fermi liquids as in 45]). 
We leave this for future investigations. 
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